Analytical advances alleviate model misspecification in non-Brownian multivariate comparative methods.

Adams and Collyer argue that contemporary multivariate (Gaussian) phylogenetic comparative methods are prone to favouring more complex models of evolution and sometimes rotation invariance can be an issue. Here we dissect the concept of rotation invariance and point out that, depending on the understanding, this can be an issue with any method that relies on numerical instead of analytical estimation approaches. We relate this to the ongoing discussion concerning phylogenetic principal component analysis. Contrary to what Adams and Collyer found, we do not observe a bias against the simpler Brownian motion process in simulations when we use the new, improved, likelihood evaluation algorithm employed by mvSLOUCH, which allows for studying much larger phylogenies and more complex model setups.

Playing Peekaboo with a Master Manipulator: Metagenetic Detection and Phylogenetic Analysis of Wolbachia Supergroups in Freshwater Invertebrates.

The infamous "master manipulators"-intracellular bacteria of the genus Wolbachia-infect a broad range of phylogenetically diverse invertebrate hosts in terrestrial ecosystems. Wolbachia has an important impact on the ecology and evolution of their host with documented effects including induced parthenogenesis, male killing, feminization, and cytoplasmic incompatibility. Nonetheless, data on Wolbachia infections in non-terrestrial invertebrates are scarce. Sampling bias and methodological limitations are some of the reasons limiting the detection of these bacteria in aquatic organisms. In this study, we present a new metagenetic method for detecting the co-occurrence of different Wolbachia strains in freshwater invertebrates host species, i.e., freshwater Arthropoda (Crustacea), Mollusca (Bivalvia), and water bears (Tardigrada) by applying NGS primers designed by us and a Python script that allows the identification of Wolbachia target sequences from the microbiome communities. We also compare the results obtained using the commonly applied NGS primers and the Sanger sequencing approach. Finally, we describe three supergroups of Wolbachia: (i) a new supergroup V identified in Crustacea and Bivalvia hosts; (ii) supergroup A identified in Crustacea, Bivalvia, and Eutardigrada hosts, and (iii) supergroup E infection in the Crustacea host microbiome community.

Model Selection Performance in Phylogenetic Comparative Methods Under Multivariate Ornstein-Uhlenbeck Models of Trait Evolution.

The advent of fast computational algorithms for phylogenetic comparative methods allows for considering multiple hypotheses concerning the co-adaptation of traits and also for studying if it is possible to distinguish between such models based on contemporary species measurements. Here we demonstrate how one can perform a study with multiple competing hypotheses using mvSLOUCH by analyzing two data sets, one concerning feeding styles and oral morphology in ungulates, and the other concerning fruit evolution in Ferula (Apiaceae). We also perform simulations to determine if it is possible to distinguish between various adaptive hypotheses. We find that Akaike's information criterion corrected for small sample size has the ability to distinguish between most pairs of considered models. However, in some cases there seems to be bias towards Brownian motion or simpler Ornstein-Uhlenbeck models. We also find that measurement error and forcing the sign of the diagonal of the drift matrix for an Ornstein-Uhlenbeck process influences identifiability capabilities. It is a cliché that some models, despite being imperfect, are more useful than others. Nonetheless, having a much larger repertoire of models will surely lead to a better understanding of the natural world, as it will allow for dissecting in what ways they are wrong. [Adaptation; AICc; model selection; multivariate Ornstein-Uhlenbeck process; multivariate phylogenetic comparative methods; mvSLOUCH.].

Controversies around the statistical presentation of data on mRNA-COVID 19 vaccine safety in pregnant women.

The work entitled "Preliminary Findings of mRNA Covid-19 Vaccine Safety in Pregnant Persons" published on April 21, 2021, in The New England Journal of Medicine, presented data collected from American surveillance systems and registries. However, problems with an unanimous interpretation of those results appeared in the public debate and citing articles. Some stated that the risk of miscarriage in vaccinated women was similar to historical values reported before the vaccines' approval. The others stated that risk was highly above-normative in women vaccinated during the first and second trimesters. We found several problems with the statistical treatment/interpretation of the originally presented values: a substantial percentage (up to 95.6%) of missing data, an incorrect denominator used for risk estimation, and too short follow-up that disabled the evaluation of the study's endpoint in numerous participants. Eventually, the Authors published a corrigendum on September 8, 2021, and pointed to updated data. Herein, we explain the statistical controversies raised by the original presentation and stress that analyzing the trade-off between knowledge and confusion brought by the release of incomplete results of such a high social interest, should aid in solving the dilemma of whether to publish preliminary data or none.

Simple SIR models with Markovian control.

We consider a random dynamical system, where the deterministic dynamics are driven by a finite-state space Markov chain. We provide a comprehensive introduction to the required mathematical apparatus and then turn to a special focus on the susceptible-infected-recovered epidemiological model with random steering. Through simulations we visualize the behaviour of the system and the effect of the high-frequency limit of the driving Markov chain. We formulate some questions and conjectures of a purely theoretical nature.

Squaring within the Colless index yields a better balance index.

The Colless index for bifurcating phylogenetic trees, introduced by Colless (1982), is defined as the sum, over all internal nodes v of the tree, of the absolute value of the difference of the sizes of the clades defined by the children of v. It is one of the most popular phylogenetic balance indices, because, in addition to measuring the balance of a tree in a very simple and intuitive way, it turns out to be one of the most powerful and discriminating phylogenetic shape indices. But it has some drawbacks. On the one hand, although its minimum value is reached at the so-called maximally balanced trees, it is almost always reached also at trees that are not maximally balanced. On the other hand, its definition as a sum of absolute values of differences makes it difficult to study analytically its distribution under probabilistic models of bifurcating phylogenetic trees. In this paper we show that if we replace in its definition the absolute values of the differences of clade sizes by the squares of these differences, all these drawbacks are overcome and the resulting index is still more powerful and discriminating than the original Colless index.

Are official confirmed cases and fatalities counts good enough to study the COVID-19 pandemic dynamics? A critical assessment through the case of Italy.

As the COVID-19 outbreak is developing the two most frequently reported statistics seem to be the raw confirmed case and case fatalities counts. Focusing on Italy, one of the hardest hit countries, we look at how these two values could be put in perspective to reflect the dynamics of the virus spread. In particular, we find that merely considering the confirmed case counts would be very misleading. The number of daily tests grows, while the daily fraction of confirmed cases to total tests has a change point. It (depending on region) generally increases with strong fluctuations till (around, depending on region) 15-22 March and then decreases linearly after. Combined with the increasing trend of daily performed tests, the raw confirmed case counts are not representative of the situation and are confounded with the sampling effort. This we observe when regressing on time the logged fraction of positive tests and for comparison the logged raw confirmed count. Hence, calibrating model parameters for this virus's dynamics should not be done based only on confirmed case counts (without rescaling by the number of tests), but take also fatalities and hospitalization count under consideration as variables not prone to be distorted by testing efforts. Furthermore, reporting statistics on the national level does not say much about the dynamics of the disease, which are taking place at the regional level. These findings are based on the official data of total death counts up to 15 April 2020 released by ISTAT and up to 10 May 2020 for the number of cases. In this work, we do not fit models but we rather investigate whether this task is possible at all. This work also informs about a new tool to collect and harmonize official statistics coming from different sources in the form of a package for the R statistical environment and presents the "COVID-19 Data Hub."

Fast likelihood calculation for multivariate Gaussian phylogenetic models with shifts.

Phylogenetic comparative methods (PCMs) have been used to study the evolution of quantitative traits in various groups of organisms, ranging from micro-organisms to animal and plant species. A common approach has been to assume a Gaussian phylogenetic model for the trait evolution along the tree, such as a branching Brownian motion (BM) or an Ornstein-Uhlenbeck (OU) process. Then, the parameters of the process have been inferred based on a given tree and trait data for the sampled species. At the heart of this inference lie multiple calculations of the model likelihood, that is, the probability density of the observed trait data, conditional on the model parameters and the tree. With the increasing availability of big phylogenetic trees, spanning hundreds to several thousand sampled species, this approach is facing a two-fold challenge. First, the assumption of a single Gaussian process governing the entire tree is not adequate in the presence of heterogeneous evolutionary forces acting in different parts of the tree. Second, big trees present a computational challenge, due to the time and memory complexity of the model likelihood calculation. Here, we explore a sub-family, denoted GLInv, of the Gaussian phylogenetic models, with the transition density exhibiting the properties that the expectation depends Linearly on the ancestral trait value and the variance is Invariant with respect to the ancestral value. We show that GLInv contains the vast majority of Gaussian models currently used in PCMs, while supporting an efficient (linear in the number of nodes) algorithm for the likelihood calculation. The algorithm supports scenarios with missing data, as well as different types of trees, including trees with polytomies and non-ultrametric trees. To account for the heterogeneity in the evolutionary forces, the algorithm supports models with "shifts" occurring at specific points in the tree. Such shifts can include changes in some or all parameters, as well as the type of the model, provided that the model remains within the GLInv family. This contrasts with most of the current implementations where, due to slow likelihood calculation, the shifts are restricted to specific parameters in a single type of model, such as the long-term selection optima of an OU process, assuming that all of its other parameters, such as evolutionary rate and selection strength, are global for the entire tree. We provide an implementation of this likelihood calculation algorithm in an accompanying R-package called PCMBase. The package has been designed as a generic library that can be integrated with existing or novel maximum likelihood or Bayesian inference tools.

Automatic generation of evolutionary hypotheses using mixed Gaussian phylogenetic models.

Phylogenetic comparative methods are widely used to understand and quantify the evolution of phenotypic traits, based on phylogenetic trees and trait measurements of extant species. Such analyses depend crucially on the underlying model. Gaussian phylogenetic models like Brownian motion and Ornstein-Uhlenbeck processes are the workhorses of modeling continuous-trait evolution. However, these models fit poorly to big trees, because they neglect the heterogeneity of the evolutionary process in different lineages of the tree. Previous works have addressed this issue by introducing shifts in the evolutionary model occurring at inferred points in the tree. However, for computational reasons, in all current implementations, these shifts are "intramodel," meaning that they allow jumps in 1 or 2 model parameters, keeping all other parameters "global" for the entire tree. There is no biological reason to restrict a shift to a single model parameter or, even, to a single type of model. Mixed Gaussian phylogenetic models (MGPMs) incorporate the idea of jointly inferring different types of Gaussian models associated with different parts of the tree. Here, we propose an approximate maximum-likelihood method for fitting MGPMs to comparative data comprising possibly incomplete measurements for several traits from extant and extinct phylogenetically linked species. We applied the method to the largest published tree of mammal species with body- and brain-mass measurements, showing strong statistical support for an MGPM with 12 distinct evolutionary regimes. Based on this result, we state a hypothesis for the evolution of the brain-body-mass allometry over the past 160 million y.

PSICA: Decision trees for probabilistic subgroup identification with categorical treatments.

Personalized medicine aims at identifying best treatments for a patient with given characteristics. It has been shown in the literature that these methods can lead to great improvements in medicine compared to traditional methods prescribing the same treatment to all patients. Subgroup identification is a branch of personalized medicine, which aims at finding subgroups of the patients with similar characteristics for which some of the investigated treatments have a better effect than the other treatments. A number of approaches based on decision trees have been proposed to identify such subgroups, but most of them focus on two-arm trials (control/treatment) while a few methods consider quantitative treatments (defined by the dose). However, no subgroup identification method exists that can predict the best treatments in a scenario with a categorical set of treatments. We propose a novel method for subgroup identification in categorical treatment scenarios. This method outputs a decision tree showing the probabilities of a given treatment being the best for a given group of patients as well as labels showing the possible best treatments. The method is implemented in an R package psica available on CRAN. In addition to a simulation study, we present an analysis of a community-based nutrition intervention trial that justifies the validity of our method.

Embracing heterogeneity: coalescing the Tree of Life and the future of phylogenomics.

Building the Tree of Life (ToL) is a major challenge of modern biology, requiring advances in cyberinfrastructure, data collection, theory, and more. Here, we argue that phylogenomics stands to benefit by embracing the many heterogeneous genomic signals emerging from the first decade of large-scale phylogenetic analysis spawned by high-throughput sequencing (HTS). Such signals include those most commonly encountered in phylogenomic datasets, such as incomplete lineage sorting, but also those reticulate processes emerging with greater frequency, such as recombination and introgression. Here we focus specifically on how phylogenetic methods can accommodate the heterogeneity incurred by such population genetic processes; we do not discuss phylogenetic methods that ignore such processes, such as concatenation or supermatrix approaches or supertrees. We suggest that methods of data acquisition and the types of markers used in phylogenomics will remain restricted until a posteriori methods of marker choice are made possible with routine whole-genome sequencing of taxa of interest. We discuss limitations and potential extensions of a model supporting innovation in phylogenomics today, the multispecies coalescent model (MSC). Macroevolutionary models that use phylogenies, such as character mapping, often ignore the heterogeneity on which building phylogenies increasingly rely and suggest that assimilating such heterogeneity is an important goal moving forward. Finally, we argue that an integrative cyberinfrastructure linking all steps of the process of building the ToL, from specimen acquisition in the field to publication and tracking of phylogenomic data, as well as a culture that values contributors at each step, are essential for progress.

Multi-omic analysis of signalling factors in inflammatory comorbidities.

BACKGROUND: Inflammation is a core element of many different, systemic and chronic diseases that usually involve an important autoimmune component. The clinical phase of inflammatory diseases is often the culmination of a long series of pathologic events that started years before. The systemic characteristics and related mechanisms could be investigated through the multi-omic comparative analysis of many inflammatory diseases. Therefore, it is important to use molecular data to study the genesis of the diseases. Here we propose a new methodology to study the relationships between inflammatory diseases and signalling molecules whose dysregulation at molecular levels could lead to systemic pathological events observed in inflammatory diseases. RESULTS: We first perform an exploratory analysis of gene expression data of a number of diseases that involve a strong inflammatory component. The comparison of gene expression between disease and healthy samples reveals the importance of members of gene families coding for signalling factors. Next, we focus on interested signalling gene families and a subset of inflammation related diseases with multi-omic features including both gene expression and DNA methylation. We introduce a phylogenetic-based multi-omic method to study the relationships between multi-omic features of inflammation related diseases by integrating gene expression, DNA methylation through sequence based phylogeny of the signalling gene families. The models of adaptations between gene expression and DNA methylation can be inferred from pre-estimated evolutionary relationship of a gene family. Members of the gene family whose expression or methylation levels significantly deviate from the model are considered as the potential disease associated genes. CONCLUSIONS: Applying the methodology to four gene families (the chemokine receptor family, the TNF receptor family, the TGF- ß gene family, the IL-17 gene family) in nine inflammation related diseases, we identify disease associated genes which exhibit significant dysregulation in gene expression or DNA methylation in the inflammation related diseases, which provides clues for functional associations between the diseases.

Exact and approximate limit behaviour of the Yule tree's cophenetic index.

In this work we study the limit distribution of an appropriately normalized cophenetic index of the pure-birth tree conditioned on n contemporary tips. We show that this normalized phylogenetic balance index is a submartingale that converges almost surely and in L2. We link our work with studies on trees without branch lengths and show that in this case the limit distribution is a contraction-type distribution, similar to the Quicksort limit distribution. In the continuous branch case we suggest approximations to the limit distribution. We propose heuristic methods of simulating from these distributions and it may be observed that these algorithms result in reasonable tails. Therefore, we propose a way based on the quantiles of the derived distributions for hypothesis testing, whether an observed phylogenetic tree is consistent with the pure-birth process. Simulating a sample by the proposed heuristics is rapid, while exact simulation (simulating the tree and then calculating the index) is a time-consuming procedure. We conduct a power study to investigate how well the cophenetic indices detect deviations from the Yule tree and apply the methodology to empirical phylogenies.

Predicting pathogenicity behavior in Escherichia coli population through a state dependent model and TRS profiling.

The Binary State Speciation and Extinction (BiSSE) model is a branching process based model that allows the diversification rates to be controlled by a binary trait. We develop a general approach, based on the BiSSE model, for predicting pathogenicity in bacterial populations from microsatellites profiling data. A comprehensive approach for predicting pathogenicity in E. coli populations is proposed using the state-dependent branching process model combined with microsatellites TRS-PCR profiling. Additionally, we have evaluated the possibility of using the BiSSE model for estimating parameters from genetic data. We analyzed a real dataset (from 251 E. coli strains) and confirmed previous biological observations demonstrating a prevalence of some virulence traits in specific bacterial sub-groups. The method may be used to predict pathogenicity of other bacterial taxa.

Using the Ornstein-Uhlenbeck process to model the evolution of interacting populations.

The Ornstein-Uhlenbeck (OU) process plays a major role in the analysis of the evolution of phenotypic traits along phylogenies. The standard OU process includes random perturbations and stabilizing selection and assumes that species evolve independently. However, evolving species may interact through various ecological process and also exchange genes especially in plants. This is particularly true if we want to study phenotypic evolution among diverging populations within species. In this work we present a straightforward statistical approach with analytical solutions that allows for the inclusion of adaptation and migration in a common phylogenetic framework, which can also be useful for studying local adaptation among populations within the same species. We furthermore present a detailed simulation study that clearly indicates the adverse effects of ignoring migration. Similarity between species due to migration could be misinterpreted as very strong convergent evolution without proper correction for these additional dependencies. Finally, we show that our model can be interpreted in terms of ecological interactions between species, providing a general framework for the evolution of traits between "interacting" species or populations.

Phylogenetic effective sample size.

In this paper I address the question-how large is a phylogenetic sample? I propose a definition of a phylogenetic effective sample size for Brownian motion and Ornstein-Uhlenbeck processes-the regression effective sample size. I discuss how mutual information can be used to define an effective sample size in the non-normal process case and compare these two definitions to an already present concept of effective sample size (the mean effective sample size). Through a simulation study I find that the AICc is robust if one corrects for the number of species or effective number of species. Lastly I discuss how the concept of the phylogenetic effective sample size can be useful for biodiversity quantification, identification of interesting clades and deciding on the importance of phylogenetic correlations.

A consistent estimator of the evolutionary rate.

We consider a branching particle system where particles reproduce according to the pure birth Yule process with the birth rate λ, conditioned on the observed number of particles to be equal to n. Particles are assumed to move independently on the real line according to the Brownian motion with the local variance σ(2). In this paper we treat n particles as a sample of related species. The spatial Brownian motion of a particle describes the development of a trait value of interest (e.g. log-body-size). We propose an unbiased estimator Rn(2) of the evolutionary rate ρ(2)=σ(2)/λ. The estimator Rn(2) is proportional to the sample variance Sn(2) computed from n trait values. We find an approximate formula for the standard error of Rn(2) based on a neat asymptotic relation for the variance of Sn(2).

Quantifying the effects of anagenetic and cladogenetic evolution.

An ongoing debate in evolutionary biology is whether phenotypic change occurs predominantly around the time of speciation or whether it instead accumulates gradually over time. In this work I propose a general framework incorporating both types of change, quantify the effects of speciational change via the correlation between species and attribute the proportion of change to each type. I discuss results of parameter estimation of Hominoid body size in this light. I derive mathematical formulae related to this problem, the probability generating functions of the number of speciation events along a randomly drawn lineage and from the most recent common ancestor of two randomly chosen tip species for a conditioned Yule tree. Additionally I obtain in closed form the variance of the distance from the root to the most recent common ancestor of two randomly chosen tip species.

Time to a single hybridization event in a group of species with unknown ancestral history.

We consider a stochastic process for the generation of species which combines a Yule process with a simple model for hybridization between pairs of co-existent species. We assume that the origin of the process, when there was one species, occurred at an unknown time in the past, and we condition the process on producing n species via the Yule process and a single hybridization event. We prove results about the distribution of the time of the hybridization event. In particular we calculate a formula for all moments and show that under various conditions, the distribution tends to an exponential with rate twice that of the birth rate for the Yule process.

A phylogenetic comparative method for studying multivariate adaptation.

Phylogenetic comparative methods have been limited in the way they model adaptation. Although some progress has been made, there are still no methods that can fully account for coadaptation between traits. Based on Ornstein-Uhlenbeck (OU) models of adaptive evolution, we present a method, with R implementation, in which multiple traits evolve both in response to each other and, as in previous OU models, to fixed or randomly evolving predictor variables. We present the interpretation of the model parameters in terms of evolutionary and optimal regressions enabling the study of allometric and adaptive relationships between traits. To illustrate the method we reanalyze a data set of antler and body-size evolution in deer (Cervidae).